function results = runGARCHmodel(y)
cmaesOn     = 1;

params      = zeros(3,1);
params(1,1) = 0.001;    %level of sigma^2
params(2,1) = 0.50;     %persistenc in sigma^2
params(3,1) = 0.01;    %ARCH coefficient in sigma^2
boundRoot   = 0.9999;
options     = optimset('Display','off','MaxIter',1e5,'MaxFunEvals',1e5,'TolFun',1e-12,'TolX',1e-12);   
[paramsNM,fvalNM]     = fminsearch(@likehood,params,options,y,boundRoot);

if cmaesOn == 1
    sigma              = 1;              %The step size
    insigma            = [0.01; 0.1; 0.01];
    opts.SigmaMax      = 5;              %The maximal value for sigma
    opts.MaxIter       = 1e6;            %The maximum number of iterations
    opts.PopSize       = 50;             %The population size
    opts.VerboseModulo = 100000;             %Display results after every 10'th iteration
    opts.TolFun        = 1e-6;           %Function tolerance
    opts.TolX          = 1e-6;           %Tolerance in the parameters
    opts.Plotting      = 'off';          %Dislpay plotting or not
    opts.Saving        = 'off';           %Saving results
    opts.Resume        = 'off';
    opts.Display       = 'off';
    warning('off')
    [paramsCMAES,fvalCMAES] = cmaes_dsgeDisplay(@likehood,params,sigma,insigma,opts,y,boundRoot);
end
if fvalCMAES < fvalNM
    paramsOpt = paramsCMAES;
else
    paramsOpt = paramsNM;
end

% Computing standard errors
numParams = length(params);
score     = zeros(length(y),numParams);
epsValue  = 1D-7;
for i=1:numParams
    params_peps = paramsOpt;
    params_peps(i,1) = params_peps(i,1) + epsValue;
    [~,~,logLike_peps]  = likehood(params_peps,y,1);
    
    params_meps = paramsOpt;
    params_meps(i,1) = params_meps(i,1) - epsValue;
    [~,~,logLike_meps]  = likehood(params_meps,y,1);

    score(:,i) = (logLike_peps-logLike_meps)/(2*epsValue);
end
A = zeros(numParams,numParams);
for k=1:length(y)
   A = A + score(k,:)'*score(k,:);
end
AVar = inv(A);

% Computing conditional vol. 
[sumLogL,sigma2Hat]  = likehood(paramsOpt,y,boundRoot);

% Computing an estimate of the variability in conditional vol, provided the
% sample is not too long (in which case it is simulated data)
rng(1,'twister')
numSim = 5000;
biasConVol2 = zeros(10,1);
if sum(~isnan(AVar(:))) && length(y)<500
    AVarChol  = chol(AVar,'lower');
    sigmaHatSim = zeros(length(y),numSim);
    sumConVol2  = zeros(10,numSim);
    s           = 1;
    while s <=numSim
        paramsMC = paramsOpt + AVarChol*randn(numParams,1);
        [~,sigma2HatMC]  = likehood(paramsMC,y,1);
        if isnan(sigma2HatMC) 
            %s
        else
            sigmaHatSim(:,s) = sqrt(sigma2HatMC);
            sumConVol2(1,s) = sigmaHatSim(:,s)'*sigmaHatSim(:,s);
            for j=1:9
                sumConVol2(1+j,s) = sigmaHatSim(1:end-j,s)'*sigmaHatSim(1+j:end,s);
            end
            s = s+1;
        end
    end
    biasConVol2(1,1) = mean(sumConVol2(1,:),2) - sqrt(sigma2Hat)'*sqrt(sigma2Hat);
    for j=1:9
        biasConVol2(1+j,1) = mean(sumConVol2(1+j,:),2)...
            - sqrt(sigma2Hat(1:end-j,1))'*sqrt(sigma2Hat(1+j:end,1));
    end
else
    sigmaHatSim = NaN(length(y),numSim);
end
    
% Saving output
results.params       = paramsOpt;
results.paramsSE     = sqrt(diag(AVar));
results.params_tstat = paramsOpt./sqrt(diag(AVar));
results.conVolSim    = sigmaHatSim;
results.conVol       = sqrt(sigma2Hat);
results.biasConVol2  = biasConVol2;
results.sumLogL      = -sumLogL;
results.AVar         = AVar;
warning('on')


    function [sumLogL,sigma2,logLike] = likehood(params,y,boundRoot)
       % test for stationarity
       if abs(params(2,1)+params(3,1)) >=boundRoot || params(1,1) < 0  ...
               || params(2,1) < 0  || params(3,1) < 0 %|| params(4,1) < 0 
           sumLogL = NaN;
           sigma2  = NaN;
           logLike = NaN;
           return;
       end
       T           = size(y,1);
       sigma2      = zeros(T,1);
       residuals   = y(1:T,1);
       %sigma2(1,1) = params(4,1);
       sigma2(1,1) = params(1,1)/(1-params(2,1)-params(3,1));
       logLike     = zeros(T,1);
       for t=1:T
           logLike(t,1)  =  - 0.5*log(sigma2(t,1))-0.5*residuals(t,1)^2/sigma2(t,1);
           if t<T
               sigma2(t+1,1) = params(1,1) + params(2,1)*sigma2(t,1)+params(3,1)*residuals(t,1)^2;
           end
       end
       % We use a minimization routine
       sumLogL = -sum(logLike);
    end

end

